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We analyze the anisotropy of turbulence in an electrically conducting fluid in the 
presence of a uniform magnetic field, for low magnetic Reynolds number, using the 
quasi-static approximation. In the linear limit, the kinetic energy of velocity compo- 
nents normal to the magnetic field decays faster than the kinetic energy of compo- 
nent along the magnetic field [Moffatt, JFM 28, 1967]. However, numerous numerical 
studies predict a different behavior, wherein the final state is characterized by dom- 
inant horizontal energy. We investigate the corresponding nonlinear phenomenon 
using Direct Numerical Simulations. The initial temporal evolution of the decaying 
flow indicates that the turbulence is very similar to the so-called "two-and-a-half- 
dimensional" flow [Montgomery & Turner, Phys. Fluids 25(2), 1982] and we offer an 
explanation for the dominance of horizontal kinetic energy. 
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I. INTRODUCTION 



In most geophysical and astrophysical flows, turbulence is affected by forces that distort 
significantly some of its scales in an anisotropic manner, such as the Lorentz force aris- 
ing from the presence of an external magnetic field in a conducting fluid. This specific 
turbulent dynamics forced by an imposed magnetic field is found in liquid metals flows, 
be they of geophysical nature — the melted iron core of the earth — or of academic interest 
in the laboratorj^. More recent laboratory experiments use sodium or gallium, whereas 
liquid sodium is used in industrial configurations (such as a French fast breeder reactor 
Superphenix). 

Generally, the motion of turbulent liquid metals is governed by magnetohydrodynamics 
(MHD): an induction equation for the fluctuating magnetic field ought to be coupled to the 
Navier-Stokes equations, and the latter equations are modified in turn by the Lorentz force, 
as a feedback from the magnetic field. In the presence of an external magnetic field, such 
MHD coupling results in new dissipative terms, of Ohmic nature, and selectively damped 
waves, the Alfven wave^. In these cases (liquid metal), the magnetic diffusivity in the 
induction equation is greater than the molecular viscosity in the Navier-Stokes equations 
{i.e. the magnetic Prandtl number is small compared to one). 

As discussed in section II, if, in addition, the magnetic Reynolds number is small enough, 
the linear regime no longer admits Alfven waves solutions, and the effect of the Lorentz 
force reduces to anisotropic ohmic (or Joule) dissipation term. In this approximation, the 
induction equation is so drastically simplified that it does not require a specific study at all. 

Forgetting the possible inhomogeneities arising from the presence of boundaries or inter- 
faces in the latter flows, the fact is that homogeneous anisotropic turbulence remains far 
less studied than isotropic turbulence. In the context of low magnetic Reynolds numbers, 
the response of initially isotropic turbulence to a static magnetic field is nonetheless doc- 
umented, in pioneering theoretical work^, many numericaP^ and experimental studied. 
One of the main properties of this kind of flow is the suppression of the three-dimensional 
motion due to anisotropic linear Joule dissipation, leading to a flow without variations in 
the direction of the imposed magnetic field. In the linear and inviscid regime, this final state 
is characterized by the following scaling of the vertical fluctuating velocity component w, 
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along the magnetic field, and the horizontal ones u and f P 



(1) 



However, recent numerical simulations of low magnetic Reynolds number turbulence show 
that the horizontal kinetic energy is dominant, at least at large scale^^. In the latter 
articles, the departure from equation ([T]) may be due to forcing schemes used to maintain 
the turbulence in a quasi-steady state. Another explanation proposed by Knaepen et a/P 
is that molecular viscosity might play a role in the decay of quasi-static MHD turbulence. 
Note that the decrease of the kinetic energy along the magnetic field has also been reported 
in a model based on the quasi-static approximation, but incorporating a non-isotropic model 
for viscous dissipatiorP. In this article, we present results of Direct Numerical Simulations 
(DNS) in order to analyze this nonlinear phenomenon, which is up to now considered as a 
restoration of isotropy, but still has to be elucidatecP^. 



II. GOVERNING EQUATIONS AND NUMERICAL SIMULATIONS 

We consider initially isotropic homogeneous turbulence in an incompressible conducting 
fluid, in which -u ~ ~ w The fluid is characterized by a kinematic viscosity u, a density 
p and a magnetic diffusivity i] = (cr/io) ^ is the electrical conductivity, fiQ the magnetic 
permeability. The initial rms velocity is uq, the integral length scale Iq- The Reynolds 
number and its magnetic counterpart are Re = {uqIq)/u ^ 1 and = (uqIq)/!] ^ 1. 
The flow is submitted to a uniform vertical magnetic field Bq scaled as Alfven speed as 
Bq = B / ^J~pJm). The ratio between the eddy turnover time and the Joule time is the magnetic 
interaction number = (B^Iq) / [rjuo). Within the quasi-static approximation, which implies 
that the asymptotic limit Rm goes to zero, but is nonetheless valid for Rm < l,^the Navier- 
Stokes equations become 

— + u-Vu = — Vp + iyV^u+ ^A~^— (2) 
at p T] oz^ 

F 

where F is the rotational part of the Lorentz force, is the inverse of the Laplacian 
operator and z the vertical coordinate, along the direction of Bq. Compressible effects are 
not taken into account here, so that V • w = 0. A pseudo-spectral method is used to solve 










FIG. 1. Volume rendering of the enstrophy (w^ is about 25% of the maximum value) using 
VAPORP. The resolution is 512^ (see section |IV[ ) . (a) Initial isotropic condition, (b) Final state 
for = 5 (the vertical correlation length is roughly one tenth of the numerical box), (c) Vertical 
plane extracted from figure (b). (d) Horizontal plane extracted from figure (b). 



equation (|2]). The velocity field is computed in a cubic box of side 27r with periodic bound- 
ary conditions using 256'^ Fourier modes. A spherical 2/3-truncation of Fourier modes is 
used to avoid aliasing and the time scheme is third-order Adams-Bashforth. The dissipa- 
tive terms are solved implicitly. The Eulerian velocity field is initialized using an isotropic 
pre-computation of eq.(|2]) in the hydrodynamic case, i.e. Bq = 0. The initial incompressible 
velocity field is a random superposition of Fourier modes distributed with a narrow-band 
kinetic energy spectrum E{k,t = 0) ~ A;^ exp(— 2(fc/fcj)^) peaked at fcj. Due to the im- 
posed magnetic field, the vertical velocity correlation lengths quickly increase. We therefore 
avoid interferences with the periodic boundary conditions by choosing ki larger than for a 
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hydrodynamic simulation. At the small-scale end of the spectrum, the minimum value of 
kraaJ-r) IS 1-2 foi all oui computatiouJI^l^ with 1^^ the Kolmogorov length scale. At the end 
of the pre- computation stage, the rms velocity is Uq = 0.45 and the integral scale Iq = 0.16 
yielding Re ^ 100. In the context of quasi-static MHD turbulence, it is not possible, at 
least using DNS, to perform high Reynolds number simulations without artificial effect of 
the periodic boundary conditions. We adopt here an intermediate configuration with mod- 
erate value of the Reynolds number and with no noticeable effect of the boundary conditions 



(see discussion on section IIIC). Note that spectral closures such as EDQNM may propose 
an interesting alternative, reaching high Reynolds numbers and forgetting about possible 
confinement effects from boundary condition^. 

The corresponding turbulent flow field is used as initial state for three different MHD 
simulations. In all of them Rm — 0.1 (hence ?7 ~ 1), so that the quasi-static approximation 
is justifiecP. Three different amplitudes of the imposed magnetic field are chosen, that 
correspond to three values of the interaction parameter: N = 1, 3 and 5. For each case, we 
perform additional "linear" simulations by neglecting the nonlinear advective term in eq.(|2]). 
The simulations are freely decaying, since forcing would hinder the natural development of 
anisotropy, which is the focus of this work. 



III. RESULTS FOR THE QUASI-STATIC CASE 
A. Kinetic energy and Reynolds stress tensor 

Figure |2] presents the ratio between vertical and horizontal energy for different interaction 
parameters. The linear state characterized by eq. ([T| is indeed observed when nonlinear 
interactions are removed. In the nonlinear simulations, the horizontal kinetic energy decays 
slower than the vertical one, which is consistent with previous observation^^. 

In order to identify the origin of this phenomenon, we study the evolution of the Reynolds 
stress tensor Rij = {ui{x)uj{x)) and its anisotropic contents bij = Rij/{2IC) — Sij/3 where 
K, is the total kinetic energy and 6ij is the Kronecker symbol. Considering the axisymmetry 
of the flow about the axis of Bq, only one of the diagonal terms is relevant!^. Figure [sj^a) 
plots the time-dependent 633 (t) for the three values of A^. The index 3 stands for the vertical 
direction. The initial value 633 (t = 0) = is characteristic of isotropic turbulence. In the 
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FIG. 2. Ratio between horizontal (perpendicular to the imposed magnetic field) and vertical (along 
Bq) kinetic energy versus time. Continuous lines correspond to nonlinear simulations whereas dot- 
ted lines correspond to linear simulations (i.e. the non-linear advective term in eq.Q is removed). 

linear regime (symbols in figure [S]^ a)), 633 quickly recovers 633 ~ 1/6 corresponding to the 
scaling in eq. ([T]). After a short initial growth — the larger N, the shorter this transient 
stage — , the nonlinear simulations exhibit a decay of 633. This indicates that the vertical 
kinetic energy decays faster than the total kinetic energy, in agreement with figure [2j This 
observation alone could lead to the incorrect conclusion that the quadratic nonlinearity 
in eq.(|2]) tends to restore 3D isotropy. A closer inspection of 633 provides a significantly 
different viewpoint, when considering the decomposition of 633 = 633 + bf^ in directional and 
polarization anisotropy contribution^^^*^, following: 



bl-l-[( e{k) - ^) sin^^d^fc 
6f3 = ^ / Z(k) sin^ e&'k 



(3) 
(4) 
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FIG. 3. Anisotropy tensor of the Reynolds stresses versus time, (a) 633 (t). (b) bl^{t) and b§^{t). 

where 6 is the polar angle between the wave vector k and the axis of symmetry, E{k) is 
the spherically-averaged kinetic energy spectrum, Z{k) is the polarization spectrum (see 
below). This decomposition is easier described in polar-spherical coordinates (see figure |4]) 
with unit toroidal vector e*^^) and poloidal vector e^^\ Due to incompressibility u(k) - k = 0, 
u = u^^^e'^^'^ + ■u^^-'e*^^-': each Fourier mode has a toroidal contribution iS^^ and a poloidal 
one M*^^^ (we drop the explicit k dependence in this paragraph). Assuming horizontal plane 
mirror symmetry, thus without helicity, the kinetic energy density is 

e = E^°' + E'"' = ^ (m(2)*^(2) + uW*uW) (5) 

and the polarization tensor density is 

Z = EP°' - E'°' . (6) 

A detailed presentation of this decomposition can be found inf^. Note that 633 is close to 
the so-called Shebalin angld^, whereas bfj is shownf^ to be exactly minus half the non- 
dimensional deviator of the "dimensionality tensor"!^. 

Figure jsj^b) shows the time evolution of 633 and bf;^. For all values of A^, 633 increases 
indicating that the kinetic energy density e{k) is not isotropically distributed among the 
different wavevector orientations. In the linear runs (symbols in figure [S]) 633 follows closely 
the nonlinear evolution (continuous lines in figure [sj. Unlike 633, a departure is observed 
for bf^: it decreases in the nonlinear simulations, but remains negligible in the linear runs. 
Thus, the decay of 633 observed in figure |3]^a) is due to a nonlinear decay of the polarization 
contribution bf^. 
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FIG. 4. Craya-Herring frame (e^^\e^'^\e^^^) in Fourier space. Fourier modes in the blue region 
contribute to E{k, 9) (eq.Q). The polar modes (i.e. ~ 0) contribute to horizontal kinetic energy, 
whereas equatorial modes (i.e. ~ '^/'^) contribute to both vertical (along e^^)) and horizontal 
(along e^-*^)) kinetic energies. 

The analysis can be even more refined when considering the scale dependence of the 
poloidal/toroidal anisotropy. In eq.([2]), one observes that the Lorentz force introduces 
an additional anisotropic dissipation, upon examining the Fourier transform of -F: F = 
— {Bq ■ hf'u/ {r]k'^) shows that the vertical modes {k/ /Bq) are strongly damped whereas 
horizontal ones (fc ± Bq) are only modified by pressure effects. This well-known phe- 
nomenon, which is linear and characterized by a Joule timescale tj = rj/B^, is responsible 
for the initial growth of 633 observed in figure |3](a). Note that tj can be used to rescale t (as 
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FIG. 5. Same as figure [sj^b) but with time normalized by the Joule timescale rj = ti/Bq. 

in figure [s]) so that all curves for 633 collapse. However, with this scaling, the evolution of 
6f3 is still very different for the various cases. As increases {i.e. the nonlinear interactions 
become more and more negligible compared to linear ohmic dissipation), the decay of bf^ is 
slower, clearly indicating that the polarization effect is triggered by a nonlinear phenomenon. 

B. Angular energy spectra 

Let us define angular energy spectra defined by 

E{k,e) = f{e) ^ (7) 

T{k,e) 

where T{k,6) denotes the torus-shaped volume defined from k — Ak/2 < |fc| < A; + Ak/2 
and 9- A9/2 < 9 < 9 + A9/2, with Ak = 1 and = vr/lO the Fourier space discretization 
steps (blue domain in figure |4]). f{9) is the geometrical weighting function such that the 
angular spectra collapse in isotropic turbulence. Torus-averaged angular spectra E{k, 9) 
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FIG. 6. Angular energy spectra (note that for the case N = 5, the results are plotted for t ~ 20tj 
whereas for the case = 1, it corresponds to t ~ 5rj). Top: Linear simulations. Bottom: Non- 
linear simulations. Left (figures (a) and (b)): = 1. Right (figures (c) and (d)) A^ = 5. On figures 
(b) and (d), the embedded figures corresponds to the same results presented in linear scale. 

have already been used in the context of rotating turbulenc^^ and are similar to the so- 
called ring decomposition^. On figures [6| we present the angular spectra for five orientations 
from the equator (6' ~ 7r/2) to the pole (6* ~ 0), for = 1 (figures [6|^a) and (b)) and = 5 
(figures [6|^c) and (d)). In addition, we distinguish the poloidal and toroidal spectra -E^"' and 
E^°^, and the linear simulations from the nonlinear ones (figures |6](a) (c) and figures [6]^b) (d) 
respectively). The figures show that almost all the energy is concentrated in the equatorial 
spectrum, as a result from the linear Joule dissipation, independently on the poloidal or 
toroidal contributions, and in both the linear and nonlinear simulations. The reasons for the 
observed differences between the latter two cases are two-fold: (a) the nonlinear downscale 
energy cascade increases the energy in the small scale range of the nonlinear runs with respect 
to the linear ones; (b ) in the meantime, a nonlinear angular transfer kicks in, explaining the 
polar energy depletion, in the same figures, for the linear runs (see e.gW). Given that polar 
modes contribute only to the horizontal kinetic energy (clearly from geometrical reasons, 
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see fig|4]) combined with a more efficient Joule dissipation in the same region, this hnear 
mechanism is responsible for the rapid decrease of horizontal energy, i.e. rapid increase of 
633, observed at short times on figjsj^a). There is however nothing in this mechanism that 
explains the decay, of bf^ in the nonlinear simulation results plotted on figjsj^b). 

The crucial difference between linear and nonlinear simulations in figures [6] is brought to 
light by the poloidal-toroidal decomposition. At the equator in the full DNS, the horizontal 
kinetic energy (also -E*™ at this specific orientation) is dominant at large scales, whereas 
the vertical kinetic energy (also E^"^) is dominant at small scales (the small plot in inset of 
figures [6|^b) and (d), in linear scale, shows more clearly the important large-scale energy gap 
between the two modes). This result explains the increase of horizontal kinetic energy with 
respect to the vertical energy, hence the decay of bf^. It is clearly of nonlinear origin since 
it disappears when nonlinear interactions are removed (see figures |6](a) and (c)). Although 
already observed in previous work^, this unpredicted dominance of horizontal energy was 
not understood. We provide here a more precise analysis since only modes such that k ± Bq 
appear to have noticeably different poloidal and toroidal energies. 



C. Correlation lengths 

We present in this paragraph some results about correlation lengths, which can be com- 
puted directly from the angular dependent spectra, and we thus check that our simulations 
are free from any artificial confinement due to periodic boundary conditions. Let us intro- 
duce the velocity correlation lengths defined by 

1 f°° 

L\. = ^ / {u,{x)u^{x + r)) dr (8) 

where Vk = r6ki is the two-point separation. In the current axisymmetric configuration, the 
most relevant anisotropy indicators involve the integral length scale with vertical separation 
but relative to either vertical or horizontal velocity componentd^ 

Ll^ = — [e(fc) + ^Z{k)] kdk (9) 
- [e(k) - ^Z(k)] kdk . (10) 

The correlation length Lf^ is characterized by the most critical growth rate, as qualitatively 
observed on visualizations [l|^b) and (d). At the end of the simulations {i.e. for t ~ 1.8), 
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FIG. 7. Evolution with time of the quantity U3L33 — lufLfi. The nonhnear results are represented 
in continuous lines and linear results in symbols. 

L^^ ^ 2, 1.6, 1 for the respective cases = 5, 3, 1. In all cases, L^^ is significantly smaller 
than the box size 27r, except for = 5 for which the largest integral scale is about one third 
of the numerical box size. However, we have seen that the results for the case = 1 and 

= 5 are qualitatively very similar (see for example figures |6|. We are therefore confident 
in the fact that our simulations are free from any spurious periodic effects. 

Moreover, it is possible to isolate the contribution due to polarization looking at the 
quantity: 



"3-^33 



2ulL\^ 



kdk 



(11) 



plotted on figure [7} This quantity has two main advantages: computing from equations 
(|9| and (10), its departure from zero is only due to the polarization Z{k); secondly, this 
quantity is accessible experimentally and is thus of particular interest. Again, one observes 
on figure [7] that polarization is negligible in the linear case whereas a significant growth is 
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observed in the nonlinear case. 

In short, the anisotropic Lorentz force is responsible for a preferential dissipation leading 
to a turbulent flow independent of the vertical direction (see figure [T]^b)); this linear effect 
concentrates the kinetic energy among horizontal modes, as attested by the growth of 633. 
Then, due to nonlinearity, the energy is anisotropically distributed among poloidal and 
toroidal components leading to a nonzero polarization and negative bf;^. 



IV. THE ANALOGY WITH 2D-3C TURBULENCE 

The quasi-static MHD turbulence is however far from being purely two-dimensional. 
Although the flow tends to be invariant in the vertical direction, the vertical kinetic energy 
is not negligible compared to the horizontal one (and is even dominant in the linear regime). 
This state is often referred to as two-dimensional, three-components (2D-3C) or "two-and- 
a-half-dimensional" turbulenc^^. We now propose a simple description of the above results 
on directional anisotropy and polarization. We first write the Navier-Stokes equations for a 
flow independent of the vertical direction z, so that d/dz = 0: 

— + ■ ^±w = uVlw (12) 
ot 

^ + u^-V±u^ = --V±p + u\/lu^ (13) 
ot p 

where u±_ = {u,v) is the horizontal velocity component and V_l the horizontal gradient. In 



both equations, the Lorentz force disappears altogether. Eq.(12) for the vertical velocity 



is that of a passive scalar, whereas eq.(13) for horizontal velocity is characteristic of 2D 
hydrodynamic turbulence. One therefore expects an inverse energy cascade for the horizontal 
velocity (or at least a strongly attenuated direct cascade) and a classical direct cascade for 
the vertical velocity. This is consistent with the results observed on figures and (d) 
since the slope of the toroidal {i.e. horizontal) component is steeper than the slope of the 
poloidal {i.e. vertical) component for horizontal modes {i.e. at the equator). Thus, we 
explain the large-scale dominance of horizontal energy as a result of 2D-like cascade for the 
horizontal velocity. The vertical component behaves like a passive scalar, and is therefore 
characterized by a direct cascade so that the vertical energy is dominant at small scales. 
Note that this effect has nothing to do with a restoration of isotropy, even if 633 tends to its 
isotropic value. 
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FIG. 8. (a) 633(4) from a 2D-3C simulation, (b) Equatorial energy spectra £'p°^'*°''(/c_l). Red/top: 
Quasi-static 3D MHD using 512'^ Fourier modes. Blue/bottom: 2D-3C simulation using 512^ 
Fourier modes (shifted down by two decades) . 



To assess the validity of this analysis, we perform two additional simulations. First, 
we compute the evolution of hydrodynamic turbulence from an initial 2D-3C state and a 
resolution of 512^. Figure |8](a) shows that 633 = 1/6 during all the simulation. As in 
quasi-static MHD turbulence, 2D-3C turbulence is characterized by a negative polarization 
indicating a dominance of toroidal kinetic energy with respect to the poloidal one. Figure 
[s^b) shows the equatorial toroidal and poloidal energy spectra E^°'^'^°^{k±) (bottom of the 
figure) that again exhibit a dominance of toroidal energy at large scales and the opposite 
at small scales, in support to our previous analysis on the final state of quasi-static MHD 
turbulence. Secondly, we would like to validate the universality of the mechanism at higher 
Reynolds number MHD turbulence, by performing a quasi-static simulation at Re ~ 300, 
using 512^ Fourier modes. We choose the case = 5 and Rm = 0.1. In this 3D simulation, 
the crossing of the toroidal and poloidal spectra (top spectra in figure [sj^b)) appears as in 
figures [6]^b) and (d) and in the previously described 2D-3C spectra. The and slopes 
are plotted for comparison with common scalings of 2D turbulence with passive scalai'^^'^, 
again showing clearly that the energy cascade is more efficient for the toroidal energy than 
for the poloidal one, resulting in a steeper slope for toroidal kinetic energy. Of course, many 
differences arise looking at figure [sj^b). This is due to the fact that quasi-static turbulence 
is not purely invariant in the vertical direction (see visualization on figjl^b)). 
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V. CONCLUSION 



In conclusion, we have found indications that the final state of quasi-static MHD tur- 
bulence, reached after a time sufficient for developping a strong anisotropy, is analogous to 
2D-3C hydrodynamic turbulence. The mechanism responsible for the transition from 3D 
isotropic turbulence to 2D-3C turbulence is subtle. It combines the anisotropic Joule dissi- 
pation, which is a linear phenomenon, and nonlinear energy transfers. If nonlinearities are 
neglected, this final state is characterized by a dominance of vertical kinetic energ How- 
ever, from the quasi-2D state arising from Joule dissipation, nonlinearity induces different 
dynamics for toroidal/horizontal velocity components and for poloidal/vertical components. 
The horizontal flow behaves like 2D turbulence whereas the vertical flow behaves like a 
passive scalar advected by 2D turbulence. Accordingly, the energy cascade (and thus the 
dissipation) is more efficient in the vertical direction, explaining the overall dominance of 
horizontal energy. We have also shown that the poloidal/toroidal decomposition and the 
distinction between directional and polarization anisotropy are fundamental to explain these 
physical phenomena. This description may help to explain the large scale anisotropy of some 
geophysical or astrophysical conducting fluid flows, with the possibility of coupling the 2D- 
3C model for the large scales to a speciflc turbulence closure at small scales, thus allowing 
to achieve very high Reynolds number simulations. 

The authors thank the computing center IDRIS of CNRS for the allocation of cpu time 
under project numbers 071433 and 022206. 
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